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ABSTRACT 

We present the results of three-dimensional radiation-hydrodynamics simulations of the formation and 
evolution of early Hn/Heiii regions around the first stars. Cooling and recollapse of the gas in the relic 
Hn region is also followed in a full cosmological context, until second-generation stars are formed. We 
£f) . first carry out ray-tracing simulations of ionizing radiation transfer from the first star. Hydrodynamics 

^ ' is directly coupled with photo-ionization heating as well as radiative and chemical cooling. The photo- 

CTN \ ionized hot gas is evacuated out of the host halo at a velocity of ~ 30 km/sec. This radiative feedback 

effect quenches further star- format ion within the halo for over tens to a hundred million years. We show 
that the thermal and chemical evolution of the photo-ionized gas in the relic Hn region is remarkably 
different from that of a neutral primordial gas. Efficient molecular hydrogen production in the recom- 
bining gas enables it to cool to ~ 100 K, where fractionation of HD/H2 occurs. The gas further cools 
by HD line cooling down to a few tens Kelvin. Interestingly, at high redshifts [z > 10), the minimum 
gas temperature is limited by that of the cosmic microwave background with Tcmb = 2.728(1 + z). The 
gas cloud goes run-away collapse when its mass is ~ 40M Q , which is significantly smaller than a typical 
clump mass of ~ 200 — 300M Q for early primordial gas clouds. We argue that massive, rather than very 
massive, primordial stars may form in the relic Hn region. Such stars might be responsible for early 
metal-enrichment of the interstellar medium from which recently discovered hyper metal-poor stars were 
C/3 , born. 

Subject headings: cosmology: theory - early universe - stars formation 
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1. INTRODUCTION 



The cosmic dark ages ended when the first sources of light turned on. These objects contributed to cosmic reionization, 
■ which observations suggest began about a few hundred million years after the Big Bang (Spergel et al. 2006; Page et al. 
2006), but the exact nature of the sources and how the process evolved are yet unknown. In the hierarchical structure 
formation model based on Cold Dark Matter (CDM), reionization is characterized by the emergence of early Hn regions 
around individual sources (stars, galaxies or quasars), followed by percolation of the ionized regions (Gnedin & Ostriker 
1997; Miralda-Escude et al. 2000; Ciardi, Ferrara & White 2003; Sokasian et al. 2003, 2004; Furlanetto et al. 2004a, 2006; 
Kuhlen & Madau 2005). The shape and extent of the early Hn regions determine the global topology of the distribution 
of neutral and ionized gas at different epochs (e.g. McQuinn et al. 2005, 2006; Zahn et al. 2005, 2006), which can be 
probed by future ground-based radio observations; turned around, the observations will provide rich information on the 
nature of the first sources of light (Madau, Meiksin & Rees 1997; Zaldarriaga et al. 2004; Furlanetto et al. 2004b; Mellcma 
et al. 2006). 

The majority of recent theoretical models indicate that a dominant contribution to the ultra-violet photons that reionized 
the Universe came from high-redshift, low-mass galaxies. These models, however, rely on the crucial assumption that 
star-formation in such low-mass galaxies is as efficient as in the observed local Universe. While such an assumption can 
reproduce the inferred Thomson optical depth to electron scattering and the neutral fraction of the intergalactic medium 
at z ~ 6, it is necessary to examine whether or not the gas within these high redshift galaxies can cool and condense 
rapidly to enable efficient star-formation. 

The standard CDM model predicts that the first cosmological objects form very early in low-mass halos in which the 
primordial gas condenses via molecular hydrogen cooling (Couchman & Rees 1986; Tegmark et al. 1997; Yoshida et al. 
2003). This mass scale is smaller than the often assumed characteristic mass for galaxy formation, for which more efficient 
atomic cooling is thought to be vital. Therefore, in the hierarchical model, feedback effects from the first generation of stars 
are expected to play a key role in setting the scene, i.e. the initial conditions, for (proto-)galaxy formation. Intriguingly, 
recent theoretical studies of the formation of the first stars (Abel, Bryan & Norman 2002; Omukai & Palla 2003; Yoshida 
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et al. 2006) consistently indicate that these objects were rather massive. It is then expected that radiation from the 
first stars would have a considerable impact on the thermal and chemical evolution of the surrounding gas cloud and 
even the intergalactic medium. In light of this, Kitayama, Yoshida, Susa & Umemura (2004, hereafter KYSU) performed 
radiation-hydrodynamics calculations of early HII regions and determined the critical mass of complete ionization of halo 
gas by the central source. KYSU showed that nearly all the gas is evacuated out of the host halo by radiation and thermal 
pressure from photo-ionization by the central massive star(s). Whalen et al. (2004) arrived at the same conclusion for 
a specific case. The long-term evolution of the relic Hii regions may critically determine the efficiency of star-formation 
in the same place, which is of considerable cosmological importance. The evolution of early Hii regions are intrinsically 
coupled with the evolution of the surrounding gas and that of the host dark matter halos, and thus are indeed a very 
complicated problem. 

There is another interesting question about the formation of the so-called second-generation stars. Unlike the somewhat 
simple and 'clean' initial conditions for the first stars, second generation stars are likely born from a gas that has been 
disturbed by earlier feedback effects. Among the most important of these effects is photo-ionization by the first stars. 
An ionized gas cools faster than it recombines, and thus many free electrons are left over (Shapiro & Kang 1987; Susa 
et al. 1998). When the temperature drops below ~ 8000K, atomic hydrogen Lyman-a cooling becomes inefficient but 
then hydrogen molecules are rapidly formed using the abundant free electrons as a catalyst. The gas cools by H 2 line 
cooling to ~ 100 K. In a primordial gas, HD molecules act as an efficient coolant at such low temperatures, enabling 
the gas to cool down to a few tens Kelvin (Flower et al. 2000). MacLow & Shull (1986) and Uehara & Inutsuka (2000) 
consider the evolution of shock-heated gas and argue that cooling by HD molecules is important under a broad range of 
conditions. The latter authors further suggested that HD cooling may lead to the formation of primordial brown dwarfs. 
Nakamura & Umemura (2002) studied the cooling of filamentary gas using an extended 14 primordial species chemistry 
that included HD and its ions. They identified a critical molecular hydrogen fraction of ~ 10~ 3 ; a gas with abundant 
hydrogen molecules can cool below 150K, and then formation, and hence cooling by HD molecules becomes important 
in the final run-away collapse phase. Interestingly, this critical fraction is close to the universal asymptotic molecular 
hydrogen fraction calculated by Susa et al. (1998) and Oh & Haiman (2002) for a cooling gas in halos with T vir > 10 4 K. 
The importance of HD cooling is further discussed in recent literature in a variety of contexts, such as in relic Hii regions 
and in supernova remnants of the first stars (Nagakura & Omukai 2005; Johnson & Bromm 2006, 2007; Yoshida 2006; 
Ripamonti 2006). Since the formation and evolution of the first Hn regions are directly linked to the formation of the 
first stars, it is necessary to perform a simulation starting from realistic initial conditions under a proper cosmological 
set-up. The relevance, and possible importance of HD cooling to early star formation can be addressed only by using such 
simulations. 

The overall sign of net radiative feedback effects from the first stars has been remaining an important, outstanding 
issue. Ricotti, Gncdin & Shull (2002) find a net positive feedback effect by ionizing radiation in initially over-dense 
regions because H 2 formation is promoted by the additional free electrons. O'Shea et al (2005) also claim an overall 
positive feedback effect, but its strength is still uncertain because their calculation does not include directly the effects of 
photo-evaporation. Oh & Haiman (2003) argue that the residual entropy of relic HII regions implies that gas collapses 
into low-mass halos with much lower central densities, making H 2 much more vulnerable to UV photo-dissociation and 
resulting in negative feedback. In particular, they predict a relation between gas entropy and the strength of the LW 
background required to quench cooling. Mesinger, Bryan & Haiman (2006) find that strong suppression of H 2 production 
persists down to the lowest redshift in their simulation when a weak soft-UV radiation background is present. This 
complicated problem clearly requires self-consistent 3D hydrodynamics and radiative transfer. 

In the present paper, we study the formation and evolution of early cosmological Hii regions using three-dimensional 
cosmological radiation-hydrodynamics simulations. Unlike previous three-dimensional calculations of Hn regions that use 
either a static density field (Alvarez et al. 2006a) or do not follow the propagation of ionization fronts (O'Shea et al. 2005), 
we couple hydrodynamics with photo-ionization heating as well as radiative and chemical cooling self-consistently. We 
first locate a primordial star-forming cloud in a large cosmological simulation. We then carry out ray-tracing calculations 
of radiation transfer. We show that nearby low-density gas clumps are destroyed by a sweeping ionization front. The 
ionized gas first escapes out of the shallow potential wells of the dark halos with a velocity much larger than the virial 
velocity, but eventually falls back when the halos have grown large enough. Consequently, there is a significant time gap 
between the formation epoch of the first generation star and that of the second one in the same comoving volume. Using 
cosmological simulations, we follow the thermal and chemical evolution of the ionized gas in relic Hn regions until the 
gas recombines, cools, and recollapses at the center of the gravitationally growing host dark halo. We show that, at high 
rcdshifts (z > 10), the minimum gas temperature of the gas cloud is limited by that of the cosmic microwave background 
with T cm b = 2.728(1 + z). Because of its low temperature, the characteristic mass of the gas clump at first run-away 
collapse is found to be ~ 4OM . Stars formed in such gas clouds likely have a smaller mass (< 40M©) than the first 
generation of stars, and hence may be progenitors of supernovae that enriched the gas from which recently discovered 
hyper metal- poor stars were born (Iwamoto et al. 2005; Johnson & Bromm 2007). Intriguingly, Tominaga et al. (2007) 
suggest that stars in this mass range trigger long-duration gamma-ray bursts with faint supernovae. 

The remainder of the paper is organized as follows. In section 2, we describe the chemistry network employed in the 
cosmological simulations. We describe our numerical techniques for radiative transfer in section 3. The results from our 
cosmological simulations are presented in section 4 and section 5. The former section describes the formation of Hn regions 
whereas the latter section discusses their evolution and the formation of second-generation stars. Finally, in section 6, we 
give concluding remarks. 
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2. PRIMORDIAL GAS CHEMISTRY 
2.1. Hydrogen and helium chemistry 

We employ the simulation code of Yoshida et al. (2003; 2006) that includes a non-equilibrium chemistry solver for 14 
species of hydrogen, helium and deuterium. We use the functional fits of Hui & Gnedin (1997) for case B recombination 
rates in the present paper. In the following, we describe in detail the relevant aspects to simulations of (relic) Hn regions. 

The rate of the charge exchange reaction 

H 2 + H+^H++H, (1) 

was recently revised by Savin et al. (2004a, b). We use this updated rate because these processes become important when 
the ionization fraction of the gas is large. Another important process is collisional dissociation of H 2 by electron impacts 



H 2 + e -> 2H + e. 



(2) 



The commonly used values are for H 2 molecules which are initially in the v = level (e.g. Shapiro & Kang 1987). Stibbc 
& Tennyson (1999) show that the LTE rate, which accounts for dissociation from higher vibrational levels, is about two 
orders of magnitude larger. Since we are interested in cooling of gas in early cosmological Hn regions, where the density 
is typically less than ~ 1cm -3 , we assume that the higher vibrational levels are not significantly populated and thus we 
use the v = rate in the present paper. 

Next, we include charge exchange between H and He 



He + +H -> Hc + H+, 



and its reverse reaction 



He + H+ 



Hc^ 



H 



(3) 



(4) 



to determine the abundance of Hc + accurately. These processes make a small but noticeable difference in He + abundance 
in a partially ionized gas at T < 5000 K. 

In a low density primordial gas, the main formation path for hydrogen molecules is via the H" channel: 



H + e -> H +hu, 
H + H- -» H 2 + e , 



(5) 
(6) 



where electrons are used as catalysts. H 2 formation via the channel and via reactions involving HeH + ions are 
dominant only at very high-redshifts (z > 200) where cosmic microwave background photons destroy by photo- 
detachment, making the H channel ineffective (see, e.g., Hirata & Padmanabhan 2006). There is considerable variation 
among published experimental data and theoretical calculations for the reaction rate of Equation (6). Glover, Savin & 
Jappsen (2006) study in detail how uncertainty in this rate affects the final H 2 abundance, and conclude that the gas 
cooling time scale is critically affected by the adopted rate. We use the fit by Galli & Palla (1998) to the calculations of 
Launay et al. (1991) for this rate, which is roughly intermediate in the range of rates studied by Glover et al. (2006). 

2.2. Deuterium Chemistry 

The chemical reactions and the rate coefficients for reactions involving deuterium are summarized in Table 1. We note 
here some differences from the previous works by Galli & Palla (1998) and Nakamura & Umemura (2002). We use the 
updated rates of Savin (2002) for the charge transfer reaction 



and its reverse reaction 

For the main HD formation path 



D + H+ -> D+ + H (reaction D2), 

D+ + H -> D + H+ (reaction D3). 
D+ + H 2 H+ + HD (reaction D7), 



(7) 

(8) 
(9) 
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Fig. 1. — The molecular cooling functions for H2 (solid), HD (dot-dashed), and IFJ (dashed). For cooling, we show both the contributions 
from e-H^~ collisions (short dashed) and H-H^" collisions (long dashed). For simplicity we assume njj = n e = 1cm -3 . 

we adopt the rate 1.6 x 10~ 9 cm 3 s^ 1 of Wang & Standi (2002) which is slightly smaller than that in Galli & Palla (1998). 
Throughout, we set the primordial deuterium abundance to the standard value of 4 x 10~ 5 . Reactions involving D~ 
are unimportant in the regime we consider, both in a neutral and in an ionized gas. Nevertheless we include them for 
completeness. 

2.3. Molecular cooling 

Radiative cooling processes owing to excitation, ionization, and recombination of atomic hydrogen and helium are well- 
determined. We use the cooling rates of Fukugita & Kawasaki (1994). We also include Compton cooling because it is 
the dominant cooling process in a diffuse ionized gas at high redshift. Below, we describe the molecular cooling processes 
that are important at low temperatures. 

We use the cooling rate of Galli & Palla (1998) for H2 line cooling, and that of Flower et al. (2000) for HD line 
cooling at low densities. Cooling by HD molecules is important only at low temperatures (T < 200K) and low densities 
(njj < 10 8 cm -3 ). At high temperatures, H2 cooling dominates because the HD abundance decreases relative to H2 (see 
section [2T4]) . The HD cooling function of Flower et al. does not include transitions between high vibrational levels but 
the contribution to the total cooling rate at the relevant densities are unimportant, as shown by Lipovka et al. (2005). 
HD cooling is effective at temperatures as low as ~ 30K. 

It is well-known that radiative cooling is limited by the cosmic microwave background (CMB) radiation at high redshifts 
because atoms and molecules act as a heating agent when T < Tcmb- The CMB temperature is given by 



Tcmb = 2.73(1 + z)K, 



(10) 



which becomes comparable to or larger than 30K at z > 10. Taking the CMB heating into account, we simply set the 
effective cooling rate as 

A = A(T gas ) - A(Tcmb). (11) 

Cooling by ionic molecules can be important in a gas with a large ionization fraction, because ionic molecules are excited 
by frequent impacts with fast electrons. We include the cooling function of H^ (Suchkov & Shchekinov 1978; Galli & 
Palla 1998). The cooling rate is given by the sum of two contributions (e-H^ collisions and H-H^ collisions), 



A c _ H + = 3.5 x 10- 27 exp T 1 „. 



A 



1.0 x 10 _28 exp 



T 
-650.0 



T n H i, 



(12) 
(13) 
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Fig. 2. — Evolution of the gas temperature (left) and molecular fractions (right) for an isobarically cooling gas. The gas is assumed to be 
fully ionized initially, with a temperature of T = 50000A". We run two cases with initial densities of njj = 1, 100cm -3 . The dashed lines in 
the left panel are for runs without HD cooling. The effect of HD cooling can be seen in the temperature and chemical evolution at t > 10 5 
years. In the right panel, we also show the fraction of H^" ionic molecules. 

where the As are energy loss rates per H^~ molecule. For the isobaric calculations we present in the next section, the H;j~ 
fraction becomes as large as ~ 10~ 7 when the temperature is ~ 5000 — 7000 K. Then the H2 fraction is ~ 10~ 4 . 

We find that the above cooling processes by H^" molecules contribute roughly equally to the total cooling rate as two 
other dominant cooling processes, hydrogen Lyman-a cooling and H2 cooling at temperatures ~ 5000 — 7000 K. We also 
consider two other ionic molecules, Hjj~ and HcH + . In addition to the fiducial chemistry network, we include reactions for 
formation and destruction of these molecules. By calculating an isobaric test case, we find that their fractions are always 
very small, and conclude that we may ignore reactions involving Hjj~ and HcH + . (See Appendix for details.) 

2.4. Isobaric cooling of an initially hot, ionized gas 

Our main objective in the present paper is not only to calculate the extent of Hn regions, but also to follow the thermal 
and chemical evolution of the ionized gas in relic Hn regions. It is intriguing that there is a (suggested) possibility of 
the formation of low mass metal-free stars in relic Hn regions (Nagakura & Omukai 2005; Vasiliev & Shchckinov 2005; 
Johnson & Bromm 2006), where enhanced molecular hydrogen production enables the gas temperature to be significantly 
lower than in a neutral primordial gas cloud. 

The evolution of an isobarically cooling gas serves as a simple, yet illustrative model of how a photo-ionized or col- 
lisionally ionized gas evolves in collapsed halos. Figure shows the temperature evolution for a parcel of gas cooling 
from T — 5 x 10 4 K. We assume that the gas is initially fully ionized and then solve the full rate equations to study 
the thermal and chemical evolution. In Figure we show the evolutionary tracks for two cases, with and without HD 
chemistry and cooling. The overall evolution except in the low temperature region appears quite similar to that in, e.g., 
Oh & Haiman (2002). We see a clear difference in the evolution at T < 200K, however. Previous calculations neglect 
the formation of and cooling by HD, and hence the gas temperature is limited to T ~ 100 K, at which point cooling by 
H2 becomes inefficient (see Fig. 1). With HD cooling, the gas further cools down to a temperature of ~ 30K, where HD 
cooling becomes inefficient. The right panel of Figure [2] shows the evolution of the H2 and HD abundances. While the 
early evolutionary tracks parallel each other, fractionation occurs at late times, when the temperature is low, and the HD 
abundance is enhanced relative to H2. The final HD abundance is found to be about one percent of the H2 abundance, 
which is a factor of 200 higher than the primordial [D/H] abundance. It is also about a factor of ten larger than that 
found in standard cosmological recombination calculations (e.g. Standi, Lepp, & Dalgarno 1998; Galli & Palla 1998). 

Fractionation of [HD/H2] occurs for various reasons. In the above example of an isobarically cooling gas, one of the 
main formation path is the reaction D + +H2 — > H + +HD, for which there is no counterpart for H2. The H2 fraction in the 
diffuse neutral gas in the early universe is ~ 10 -5 — 10~ 6 (Standi et al. 1998; Galli & Palla 1998), whereas the universal 
H2 fraction in the collapsing halo is as large as a few times 10 -3 (Susa et al. 1998; Oh & Haiman 2002). Hence, the large 
H2 fraction promotes the formation of HD molecules at low temperatures. In other words, HD is energetically favored 
because the binding energy of HD molecules is higher than that of H2 molecules. The relative abundance in equilibrium 
is 

"(HP) _n(D) /465K\ 

M^) =2 MB) exp [—) ' (14) 



and thus HD molecules are preferentially produced at temperatures much lower than 465 K (Solomon & Woolf 1973). 
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Fig. 3. — Schematic diagram showing how evaluation points are defined for a given particle distribution. The source and the target particle 
are connected by a straight line. The evaluation points (large colored circles), at which the local densities are calculated, are defined following a 
Smoothed Particle Hydrodynamics procedure (i.e. the "gather" formalism of Hernquist & Katz 1989). The spacing of two adjacent evaluation 
points is given by a fraction of the local smoothing length. For the radiative transfer calculation, we compute the ionization-front crossing 
time using the density p; and the local path length i;. 



In this section, we describe our numerical scheme for radiative transfer coupled with hydrodynamics and chemistry. 
The entire procedure is fully parallelized and coupled with the parallel smoothed particle hydrodynamics (SPH) code 
GADGET2 (Springel 2005), so that the radiative transfer calculation, as well as the gravity and hydrodynamics, can be 
performed on massively parallel architectures. We note that the version of SPH we adopt is a fully conservative formulation 
(Springel & Hernquist 2002), where the equations of motion properly account for evolution in the smoothing lengths of 
the gas particles. At the end of this section, we present a simple, yet very important, test case of photo-evaporation of a 
small gas sphere using our code. 



Our radiative transfer calculations are done in a two-step manner as follows. For a given point radiation source, we 
first compute and assign photon arrival times for all the gas particles surrounding it. The photon arrival time for a gas 
particle is calculated by integrating the ionizing photon consumption over the path from the central source to the gas 
particle. We define a path, a photon ray, by a straight line connecting the source and a target particle. We then compute 
local densities and smoothing lengths at many evaluation points on the path in an SPH fashion. Figure [3] is a schematic 
diagram showing how evaluation points are defined and the local densities are calculated. Starting from the central source, 
an ionization front is advanced on the path over the segment length U between the i-th and i + 1-th evaluation points. We 
calculate the local density pi and the corresponding smoothing length hi at the i-th point as in SPH. The position of the 
next i + 1-th point is determined by the condition that the length li is smaller than the local smoothing length hi\ i.e., 
the density variation over li must be sufficiently small. We take li = l/5hi, after performing several tests to confirm that 
halving or doubling the numerical factor does not change the result. Between two adjacent evaluation points, the I-front 
crossing time is computed by integrating the jump condition for the I-front position, 



3. RADIATIVE TRANSFER OF IONIZING PHOTONS 



3.1. Ionization front propagation 
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Fig. 4. — Characteristic time evolution of the shock radius and the I-front radius in a gas sphere with a steep density profile. The I-front 
changes from D-type to R-type at a time T = T c . The transition point -Rstrom can be evaluated from the initial Stromgren radius given by 
Equation ifTTt . 

We approximate the above integral by a summation over ray segments as 



assuming complete ionization for hydrogen and/or helium depending on the relative extents of Hn and Hem regions (see 
Section 1^1)) . We advance the I-front by repeating this procedure until we reach the target particle, or until the source 
lifetime elapses. In this way every particle is assigned its own photon arrival time t a .i. Particles having t^i < tiif e are to 
be irradiated by the radiation source at t^. 



While the above method works well when I-front propagation is not density-bounded, i.e., of R-type, there is a situation 
where it does not provide a physically correct solution for t a { (Susa 2006; Ahn & Shapiro 2007). When an I-front is of 
D-type, it can move forward only at the velocity of the foregoing hydrodynamic shock. Namely, the I-front catches up 
with the the foregoing shock as soon as the shock reduces the density inside it (see Shu 1992 and KYSU for a detailed 
description of D-type and R-type ionization fronts). For initially steep power-law density profiles, p oc r~ w with w > 1.5, 
an I-front begins as D-type, and then becomes R-type when the internal density is reduced by a foregoing hydrodynamic 
shock that erases the initial steep density profile. Through a series of numerical simulations, KYSU found that, for a given 
density profile, the time when the transition of D- to R-type occurs can be accurately estimated a priori by a Stromgren 
analysis for the initial density profile. We describe the method below. 

The Stromgren radius is obtained by equating the photon production rate of the central source and the total recombi- 
nation rate, and the result is expressed as 



where Ni on is the photon production rate of the source and n(-Rstrom) is the mean density within i?strom- The mean 
density n within the shock decreases as it propagates and sweeps outward. i?strom increases then, and eventually the 
condition i?strom = -Rshock is met. At this time, ionization is balanced by total recombination within i? s hock- This is the 
time when the I-front can take off horn, the shock and changes into R-type. Figure 2] schematically shows how the shock 
and the I-front propagate. KYSU studied the time and position of the I-front transition using one-dimensional calculations 
with full radiative transfer, and showed that equation p7|) indeed provides an accurate estimate for the transition point. 

In order to reproduce the D-to-R transition in three-dimensional calculations, we first evaluate the initial Stromgren 
radius using Equation (|17[) for the density profile around the source, and assume that the I-front moves at a constant 
shock velocity while it is of D-type. From the results of KYSU, we set the shock velocity v s hock = 25kms~ 1 . Although the 
exact value of i>shock is, in principle, dependent on the density profile and the infalling gas velocity, it is a weak function of 
the density slope over the range of gas density profiles considered in KYSU, and is about 25 — SSkms^ 1 (see also Shu et al. 
2002). Photon arrival times within R s are simply given by the time when the shock arrives. Once the 1-front transforms 
to R-type, we can use Equation (fl5| to compute photon arrival times as described in the previous section. 




(16) 



3.2. Treatment of D-type ionization front 




pc, 



(17) 
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3.3. Photo-ionizationa and heating 

The second step is the calculation of photo-ionization and photo-heating rates. Gas particles are irradiated by the 
central source after t a ^. We use an approximation that, after t a) i, the radiation intensity and spectrum at each particle's 
position are computed in the optically-thin limit. Although this assumption might seem overly simplified, it is indeed a 
good approximation for the particular problems we are interested in, where the initial gas density profile is steep. Three- 
dimensional simulations (Abel et al. 2002; Yoshida et al. 2003, 2006) as well as one-dimensional calculations (Omukai & 
Nishi 1998; Ripamonti et al. 2002) show that primordial gas clouds have steep density profiles p oc r~ n with n ~ 2 — 2.2. 
For such profiles, only the very central part is dense, where ionizing photons are effectively consumed by recombinations 
initially. This dense region is, however, quickly swept by the hydrodynamic shock and the mean recombination rate 
rapidly decreases as well. Consequently, the opacity of the interior is progressively reduced and becomes less and less 
important as time elapses. Note also that stellar radiation with 13.6eV < hv < 54.4eV is not significantly absorbed by 
HI in a Hll /Hen /Hem region for the reasons explained in section In section \WM we run a test problem and show 
that the above approximation works rather well for our particular situation. 

The photo-ionization rate for hydrogen H + hv — > H + + e~ is computed as 



d n I du^, (is) 



and the corresponding heating rate is given by 



hv 



r = n H i dQ dv^-^(hv-hv L ), (19) 

J Ju^ hv 

where hv^ is the ionization threshold energy, 13.59 eV. The photo-ionization and heating rates for helium are computed 
similarly. We adopt the ionization cross sections in Osterbrock & Ferland (2006). When helium ionization is included, as 
described in the next section, we modify the radiation intensity according to the ionization structure. 

3.4. Helium ionization 

Massive metal-free stars emit a large number of photons with hv > 24.5 eV, and those with hv > 54.4 eV as well, 
and thus can ionize helium. The relative production rate of these high energy photons (e.g., Bromm, Kudritzki & Loeb 
2001) indicates that the structure of the ionized region around a massive Population III star may be similar to that of 
present-day planetary nebulae, where a small Hem region in the immediate vicinity of the star is surrounded by a larger 
Hn/Hen region. It is important to include helium ionization because the temperature of the ionized gas is raised by 
additional photo-heating, and, consequently, the pressure gradient (the driving force of the gas flow) is increased. To 
include helium ionization, we adopt approximations suggested by Osterbrock & Ferland (2006): (1) the gas is almost fully 
ionized within the Hem region, and (2) hydrogen in the Hem region is kept ionized by recombination of Hen Lyman-a 
and Balmer photons as well as the continuum produced via Hen two-photon process. Because the number fraction ratio 
of He/H is small, the Hn and Hen regions have the same extent if the source spectrum is a hot (T ~ 10 5 K) thermal one 
such as that for massive metal- free stars (e.g., Schaerer 2002). These approximations considerably simplify the remaining 
task. 

We first compute photon- arrival times for the Hem region using Equation (|15[) . assuming the gas is fully (both hydrogen 
and helium) ionized within it. (Note that n p needs to be replaced by n- [ic ++ in equation |15|). We then compute the 
extent of the Hn region surrounding the Hem region following a similar procedure. Here, we make use of the above 
approximation (2), and compute the arrival times of hydrogen ionizing photons assuming that they are not absorbed in 
the Hem region. In general, the boundary of an Hem region may not be as sharp as we assume here, because a fraction 
of photons in the high-energy tail can penetrate further. Nevertheless, we adopt the above simplifications to obtain the 
approximate size of the Hem region. Detailed structure of helium ionization will be presented elsewhere (Kitayama et al, 
in preparation). In summary, after all of these procedures, every gas particle is assigned two photon arrival times, t^hvi 
and i a ,i,/ti/ 2 where hv\ = 13.6 eV and hv^ = 54.4 eV, respectively. The particles are irradiated by photons with hv > hv\ 
at t a .ij lVl , and by those with hv > hv2 at £ ai i,/ii> 2 ■ The photo-heating rates are calculated consistently with this. Within 
a Hn/Hen region, the radiation spectrum is cut-off above hv^ = 54.4eV, whereas within a Hii/Hem region, we use a full 
thermal spectrum for the assumed hot stellar source. 

3.5. Photo- detachment and photo-dissociation 

We include photo-detachment of H by photons with hv > 0.755 eV, and photo-dissociation of H2 by Lyman- Werner 
(LW) photons. For photo-detachment 

H-+ 7 ^H + e, (20) 

we use the cross-section of de Jong (1972) 

cr H - = 7.928 x 10 5 (^ - v th ) 3/2 v~ 3 cm 2 (21) 

for hv > hv t h — 0.755eV. Since the optical depth for these photons is generally small in the primordial intergalactic medium 
(see Reed et al. 2005 for realistic values around cosmological density peaks), we adopt the optically-thin approximation. 
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Namely, when a radiation source is turned on, the radiation intensity at 0.755eV < hv < 11.18eV is calculated in its 
optically-thin limit at every point in the simulation region. 

Photo-dissociation of hydrogen molecules by LW band photons 

H 2 + 7 -> H; ^ H + H, (22) 

is handled in the same manner. We use the dissociation rate 

few = 1.1 x 10 8 F LW s _1 (23) 

where Flw is the LW photon flux at hv — 12.87eV in units of ergs _1 cm _2 Hz~ 1 (Abel et al. 1997). This reference 
photon frequency corresponds to transitions to the v = 13 vibrational level of the Lyman state, that lies in the middle 
of the strongest transitions. To compute the LW radiation flux, we again adopt the optically-thin approximation and 
set it as Flw oc 1/i? 2 where R is the distance from the source. Note that photo-dissociation is unimportant within Hn 
regions where hydrogen molecules are destroyed mainly by collisional dissociation and charge transfer. Because of the 
small frequency range of the H2 LW band (11.18 - 13.6 eV), Equation (|23[) is a good approximation even if the radiation 
intensity varies smoothly over this range. However, the hydrogen Lyman-series absorption causes a substantial intensity 
reduction at the line frequencies, some of which lie close to a few of the strongest LW lines (Haiman, Rees & Loeb 1997; 
Ricotti et al. 2002; Yoshida et al. 2003). We ignore the effect for computational efficiency, and hence it is expected that 
photo-dissociation outside the H11 region would be slightly less effective than reported in the present paper. 

In and around dense molecular gas clouds, the effect of self-shielding needs to be taken into account. For a static 
isothermal gas, an effective shielding factor is given as a function of column density, by 



/shield = min 



10 14 cm" 2 



(24) 



where N^ 2 is the molecular hydrogen column density (Drain & Bertoldi 1996). We account for H 2 self-shielding as in 
Yoshida et al. (2003) and Glover, Savin, & Jappsen (2006). Briefly, the radiation intensity at each position in the 
simulated region is computed by assuming that it is attenuated through surrounding dense gas clouds. We define a local 
molecular hydrogen column density iVjj 2 in a consistent manner employing the SPH formalism. We line-integrate the 
molecular hydrogen number density around the z-th gas particle according to 

Nn 2 ,i = J n H2 dl, (25) 

where is the position of the i-th gas particle and r max is the length scale we choose in evaluating the column density. 
In practice, we select an arbitrary line-of-sight and sum the contributions from neighboring gas particles within r max by 
projecting an SPH spline kernel for all the particles whose volume intersects the sight-line. We repeat this procedure in 
six directions along orthogonal axes centered at the position of the i— th particle. We then take the minimum column 
density as an conservative estimate. While the shielding effect becomes significant for 7Vh 2 ^> 10 14 cm -2 , the gas remains 
nearly optically thin to LW photons even for column densities N-r 2 ~ 10 20 ~ 21 cm~ 2 if there are large velocity gradients 
and/or disordered motion (Glover & Brand 2001). When computing the above integral, we discard gas particles that 
have a large relative velocity (V Ic i 3> ^thermal) to the i-th particle, in order to account for the reduction of the shielding 
effect owing to Doppler shifting, following Glover et al. (2006). We use a constant thermal velocity ^thermal = 2 km/sec, 
by noting that an H2 rich dense gas has a low temperature of T < 500 K. We have chosen the length scale r max = 100 

physical parsec by noting that the virial radius of an early mini-halo with mass 10 6 M Q is about 100 parsec. There are 
not significantly larger gas clumps than this scale in the simulated region. The local column density estimates are easily 
computed along with other SPH variables with a small number of additional operations. 

We treat photo-dissociation of HD molecules similarly. Since some of the HD Lyman-bands lie close to those of H2 , 
the HD lines can be shielded by H2 . Barsuhn (1977) argues that the HD dissociation rate can be decreased by as much 
as a factor of three. Since HD molecules are formed in H2-rich regions, we assume that the shielding by H2 is maximally 
effective, and adopt this reduction factor for HD photo-dissociation. 

^^Here we consider the LW photons emitted directly from the central radiation source (a star) and not from the UV background radiation, 
which is considered separately in section [5]3] 
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Fig. 5. — Minihalo evaporation test problem. We plot the radial profiles of density (top-left), temperature (top-right), neutral hydrogen 
fraction (bottom-left), and radial velocity (bottom-right) at two output times, t=1.8 X 10 5 yrs and t = 2.2 X 10 6 yrs. We compare the results 
of our three-dimensional calculation (shown as grey contours) with those of the full-radiative calculation of KYSU (dashed lines). 



3.6. A test problem: Mini-halo evaporation 

We test the overall accuracy of our radiation hydrodynamics calculation using a spherical cloud problem. Expansion of 
an I-front in a spherical gas cloud offers a simple, yet illustrative case for this purpose. We adopt a spherical gas cloud 
problem studied by KYSU: a massive Population III star with M» = 200Af© is embedded at the center of a 10 6 M Q dark 
halo that is collapsing at z = 20. We set the initial gas density profile to follow a power law n = 3000(r/pc)~ 2 cm~ 3 . 
We distribute 0.2 million particles in the sphere according to this density profile. The source is chosen to have a thermal 
spectrum with effective temperature T e g = 10 5 K. For this test problem we include only hydrogen. (See KYSU for further 
details.) 

Figure [5] shows the evolution of the Hn region around the central star. We compare the radial profiles of density 
(top- left), temperature (top-right), neutral hydrogen fraction (bottom-left), and gas velocity (bottom-right) with those of 
KYSU. While there is a slight off-set in the position of the D-type I-front at t = 0.18 Myrs (note that the difference is 
only ~1 pc), the agreement between the profiles is very good both in amplitude and shape. The off-set likely owes to a 
slight difference between the assumed velocity of the D-type front and the actual shock velocity. Although we can easily 
fix the velocity of the D-type I-front so that the position of the D-type front is perfectly matched to the result of the 
ID calculation at a given epoch, we do not attempt to do so here because the difference is only transient, and the initial 
behavior of the I-front does not much affect the final results for the long-term evolution of Hn regions. In Figure [5l the 
sharp transition at r ~ 4pc (t = 0.18 Myrs) is a D-type front, which changes to R-type at t ~ 0.2 Myrs. The profiles at 
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t = 2.2 Myrs show a clear bump at r ~ 100 pc that is a shock front. The R-type ionization front is located about 1 kpc 
away from the center, so the entire region is ionized by this time, as can be seen in the neutral fraction profile. At t = 2.2 
Myrs, the radial velocity peaks at ~30 km/sec, which is also well reproduced. 

From Figure O we conclude that the accuracy of our code is quite satisfactory. The transition from D-type to R-type 
ionization front as well as hydrodynamic variables such as density and temperature are accurately reproduced in the 3-D 
calculation. Furthermore, even the neutral fraction profile is reproduced rather well. We emphasize that it is generally 
easy to obtain the correct ionized fraction within Hn regions because it is close to unity, whereas it is hard to produce a 
correct neutral fraction because it usually has a small value and is sensitive to the local density, temperature, and their 
time evolution. Hence, the level of agreement with the accurate solutions shown is reassuring. Our radiation transfer 
scheme is similar to that of Susa (2006) and the accuracy of our code to the particular problem of minihalo evaporation 
appears as good as Susa's result. We mention that Susa's code can be applied to more general problems, because it 
updates gas opacities more frequently as density structures evolve. The method described in this section is suitable for 
radiative transfer problems where there are physically separated point sources which have short lifetimes, such as those 
for early Hn regions. Further numerical implementations will be necessary for more general problems. 

4. COSMOLOGICAL SIMULATION 

We can discuss the long-term evolution of early Hn /Hem regions only if they are studied in a fully cosmological 
context. To this end, we run a cosmological simulation adopting the concordance ACDM cosmology with matter density 
f2 m = 0.26, baryon density fib = 0.04, cosmological constant fl\ = 0.7 and expansion rate at the present time Ho = 70 
km s -1 Mpc -1 . The density fluctuation amplitude is normalized by setting erg = 0.9. While some of these cosmological 
parameters are different from those determined most recently from the third- year WMAP data (Spergel et al. 2006), all of 
our results presented below are robust to slight changes in the values. Nevertheless, it is probably important to note how 
adopting a smaller value of erg would affect the epoch of the events such as first-star formation. Generally, the formation 
epoch is shifted to lower redshifts for lower values of erg. For the value from the third-year data, erg = 0.74, everything 
occurs later, about 40% in cosmological expansion parameter (e.g., Alvarez et al. 2006b). Hence, if the first stars are 
formed at z ~ 20 in the old ACDM model with erg = 0.9, then it would be at z ~ 12 — 13 in models with erg = 0.74. This 
may be important when we consider the evolution of gas in relic Hn regions. The temperature of the cosmic microwave 
background radiation, which is a function of redshift, will affect the evolution of the prestellar gas cloud collapsing in the 
relic Hn region. We discuss this issue in section [S~2l We also note that the fact that structure formation occurs late in 
low ag models may actually provide better prospects for direct observations of the first stars and early Hn regions. 

We start from a low- resolution simulation of boxsize 1.4Mpc on a side. We identify one of the most massive halos 
in the simulation volume at a redshift of 20, and then make zoomed initial conditions around this object with a much 
higher mass resolution (see Fig. [5]). In the highest resolution region, the gas particle mass is 1.82M Q and the dark matter 
particle mass is 11.8M . We use the parallel Tree N-body/SPH solver GADGET2 (Springel 2005) to evolve the system 
from z — 100 to the epoch when the first collapsed gas cloud is identified. The code includes all the relevant atomic and 
molecular processes as described in Yoshida et al. (2003, 2006). Figure [5] shows the simulated high-resolution region and 
a boundary low-resolution region. We locate a cold (T < 500 K) and dense (n > 10 3 cm -3 ) cloud whose mass exceeds 
300M Q as a site for star-formation. The critical mass is obtained from simulations of Abel et al. (2002) and Yoshida et 
al. (2006) for gas clouds at the first run-away collapse. The dark matter halo which hosts the star-forming cloud has a 
mass of 5 x 10 5 M Q . We define the halo virial radius such that the mean mass density within the radius is 200 times the 
critical density and the halo mass is the mass within the virial radius. We embed a IOOMq PopIII star at the densest 
part of the gas cloud and switch on radiation. The stellar mass is derived from the recent work of Yoshida et al. (2006) 
that follows the proto-stellar evolution in detail. For the source spectrum, we use a blackbody with effective temperature 
T c g which is determined from the stellar mass (Schaerer 2002). 

We do not trigger a supernova explosion in the simulations in the present paper for several reasons. First, the fate of 
Population III stars with mass ~ 20 — 100M Q is rather uncertain in terms of their mechanical feedback efficiency (Umcda, 
Nomoto & Nakamura 2000; Heger & Woosley 2002), while very massive stars with ~ 140 — 260M© likely cause complete 
destruction of the surrounding structure (Bromm, Yoshida & Hernquist 2003). Kitayama & Yoshida (2005) show that the 
halo destruction efficiency of supernova explosions is sensitively dependent on the explosion energy, initial density of the 
explosion site and even on reshift, because of efficient Compton cooling in the early Universe. Therefore, it is necessary to 
use outputs of simulations of Hn regions, such as those presented in the present paper, for the study of supernova feedback 
effect. We will defer a consideration of these effects to future work. 

4.1. Formation of the First HII Region 

Figure [6] shows the extent of the Hn region around the first star at t =1.1, 3 Myrs. The Hn region has a few kilo-parsec 
diameter at its maximum extent, and the ionized gas has a temperature of ~ 1 — 3 x 10 4 K (see Fig. [7]). The diameter 
is similar to what has been found in previous works on Hn regions around massive primordial stars (KYSU; Alvarez et 
al. 2006a; Abel, Wise & Bryan 2006). The Hn region has a characteristic "butterfly" shape which occurs because fast 
recombination in dense regions such as filaments prevents the I-front from propagating. Although it is hard to sec in the 
figure, a close look reveals that I-front propagation is also delayed by forming gas clumps near the main halo. However, 
since the maximum density of such small gas clumps is still ~ 1 cm -3 , they are quickly photo-evaporated within the 
central star's lifetime (Shapiro, Iliev, & Raga 2004; Ahn & Shapiro 2007; Susa & Umemura 2006). 
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Matter distribution 




Fig. 6. — Expansion of the HII region around the first star. The top-left panel shows the large-scale density distribution around the first 
object at 2 = 26. A three-level zoom-in technique was used to make the initial conditions for this simulation. Panels on the right show the 
projected gas density weighted by ionization fraction. Each panel has a side-length of about 7 proper kpc. We show the extension of the Ho 
region at t=l.l, 3 Myrs after the central star turns on, as indicated in the figure. The bottom-left panel shows the extent of the Hn region 
(green) and that of the Hem region (blue) at t = 3 Myrs. We use a different weighting and color scale for this panel. 

This is because there is only one halo in which the gas has cooled by H2 cooling at z = 26. Other surrounding halos have 
masses smaller than ~ 10 5 Af Q , which is much smaller than the critical collapse mass of ~ 5 x 10 5 M©. 

Since the gas in the main halo had initially a steep density gradient, (almost) prompt photo-ionization by the R-type 
front results in a significant pressure gradient across the H11 region, and the thermal pressure drives a strong gas wind 
outwards. The darkest round part at the center of each panel on the right in Figure [6] is the region where the outgoing 
shock is evacuating the gas out of the host halo. We show the detailed structure of the ionized region in the bottom-left 
panel. We see a smaller Hem region, of diameter ~ 1 kpc (colored in blue), surrounded by a larger Hn (and Hen ) region. 
Both hydrogen and helium in the main halo gas are nearly fully ionized, as can be seen from the extent of the Hem region. 
Within the Hem region, the gas has a temperature of ~ 3 x 10 4 K. 

Figure [7] shows the radial profiles for density, velocity, and temperature at t = and t = 3 Myrs. Thermal pressure 
quickly suppresses the steep density profile and turns it to a nearly flat one within the stellar lifetime. The ionized region 
has a high temperature of T ~ 30000K in the case with helium. For reference, we also show the temperature profile 
from a run in which we include only hydrogen. Additional heating by helium ionization is clearly seen. Throughout the 
evolution, the peak gas velocity is about 30 km/sec, which is nearly 10 times larger than the escape velocity of the host 
halo. Hence it is expected that almost all the gas will be evacuated within a few tens of million years (for a halo with 
virial radius ~ 100 pc). 
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Fig. 7. — Radial profiles of density (top), velocity (middle), and temperature (bottom) at t = and t = 3 Myrs. For the temperature 
profile, we show the results of two cases: one with helium ionization and the other with hydrogen only. (See text in section [4. 1 1 ) 
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Fig. 8. — Time-averaged ionizing photon escape fraction as a function of stellar mass. 

Because of this efficient evacuation by radiation, further star-formation in the same halo is quenched for a long time until 
the gas can cool, be re-incorporated, and condense again. We study the evolution of the gas in the relic Hn region further 
in the next section. 

An important quantity often used in simulations and theoretical models of reionization is the ionizing photon escape 
fraction. Our three-dimensional radiation hydrodynamics calculation allows to directly measure this quantity. Since the 
mass resolution of our simulations is sufficient to resolve the smallest gas clumps within early mini-halos, the effect of gas 
clumping is self-consistently included in our calculations. Therefore, we do not need to assume the so-called photon escape 
fraction nor gas clumping factor. The photon escape fraction is indeed obtained as a result of our simulations. We define 
the escape fraction as the number fraction of photons that escape out of the host halo's virial radius into the intergalactic 
medium. Following Alvarez et al. (2006a), we evaluate the escape fraction along a photon ray as 



/esc(t) = 1 



47TO!B 



P h 



n 2 (r, t) r 2 dr , 



(26) 



where r v - lr is the virial radius of the host halo. The halo has r v i r = 80 pc at z — 26. We compute the escape fraction 
averaged over the lifetime of the star. The value obtained is / esc = 0.77 for our fiducial model with M* = 100M©. Figure 
[8] show the ionizing photon escape fraction as a function of the central stellar mass. The result is consistent with the 
results of previous one-dimensional calculations assuming spherical symmetry (KYSU). We calculate a few other cases by 
varying the mass of the star. Figure [5] shows / csc as a function of the stellar mass. We find that the escape fraction as a 
function of the stellar mass can be well-fitted by 



/csc = 2.0 - exp 



1.0 



(M* - 25.0) 



0.37 



(27) 



for the particular host halo we simulate, i.e, a halo with mass 5 x 10 5 M Q at z = 26. The obtained escape fraction is 
in excellent agreement with that of Alvarez et al. (2006a); we have confirmed the result with a ~70 times better mass 
resolution. The dependence of the escape fraction on halo mass and collapse epoch, gas density profile etc. can be easily 
inferred following the description of D-to-R type transition of the I-front as discussed by KYSU. 
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Fig. 9. — Underlying dark matter density field at z = 26 when the first star is formed (left), and at z = 16 when the second-generation 
star(s) are formed (right), showing a significant evolution in ~ 100 Myrs. The main halo is marked by A. Halo B is a nearby massive halo, 
whose evolution is studied in section [5. II The descendant of Halo A is indicated in the right panel, where there is another halo (C) close to 
it. Halo C hosts the second generation star-forming gas cloud which was formed from an ionized gas. 

5. EVOLUTION OF RELIC Hn REGIONS AND THEIR SURROUNDINGS 

In this section, we study the thermal evolution and re-collapse of the gas in relic Hn regions. We first discuss possible 
star-formation in regions outside the Hn region during and after the central star's lifetime, and then show the results for 
the chemo-thermal evolution of the ionized gas. We clarify how the gas cloud evolution differs from that for primordial 
gas clouds (e.g., Yoshida et al. 2003, 2006). 

5.1. Evolution of Nearby Halos 

There are many low-mass halos in the regions surrounding the first star, where ionizing radiation does not arrive within 
the lifetime of the central star. Although the photo-dissociation region (PDR) extends much farther than the Hn region, 
the central massive star is short-lived and thus the effect of photo-dissociation is only temporal. We first study the 
impact of radiation from the first star only, ignoring the effect of a global cosmic UV background radiation field. After 
the first star dies, some of the nearby halos are able to re-form enough hydrogen molecules to meet the baryonic cooling 
criteria (Tegmark et al. 1997; Yoshida et al. 2003; Oh & Haiman 2003). Throughout the evolution, photo-dissociation 
of H2 , photo-detachment of and H2 molecule formation act in a complex manner. The Lyman- Werner photons 
continuously destroy hydrogen molecules within a large region, and photons with energy greater than 0.755eV effectively 
destroy H~(see section |3~5)) . However, the effect is weak in regions sufficiently distant from the radiation source, where 
collisional detachment processes still regulate the abundance of H~ just as in the case without radiation. During the 
lifetime of the first star, the outer PDR remains H2 poor but the H - abundance is not greatly reduced, owing to the 
difference in their respective cross-sections. The gas in distant mini-halos has a high temperature (~ T v - lr ) because of the 
absence of an effective coolant. Hence, as soon as photo-dissociating radiation is switched-off when the central star dies, 
H2 molecules are quickly reformed. It is interesting that the reformation takes place in a synchronized way, determined by 
the death of the first star. This H2 regeneration mechanism differs from the often claimed positive feedback, in which H2 
molecules are quickly formed in a (partially) ionized region behind the I-front (Ricotti, Gnedin & Shull 2002) or behind 
a M-type shock (Ahn & Shapiro 2007; Susa & Umemura 2006) 

In our simulation, there is a 'second' halo at ~ 2 kpc away, in which the gas cools and condenses after the first star dies 
(see Halo B in Fig. [5]). Interestingly, Halo B has a slightly larger mass than the main halo (Halo A) at z = 26, but the 
growth of Halo B has been so rapid that significant dynamical heating associated with progenitor mergers prevented the 
gas in Halo B from cooling, as often found in simulations of early structure formation (Yoshida et al. 2003; Reed et al. 
2005; O'Shea & Norman 2007). Hence the first star is formed earlier in Halo A. At the position of Halo B, the dissociation 
(owing to the stellar radiation from Halo A) timescale is just about 1 Myrs (from equation [23] )• The H2 fraction first 
slightly decreases, but rapid re-formation of H2 leads to efficient cooling and condensation within a few million years after 
the first star dies. Thus star-formation will take place in Halo B on a similar timescale, just as for Halo A, and the formed 
star, if massive, will produce a large Hn region. We have carried out ray-tracing calculations assuming that a 100 M Q 
star is formed in Halo B. We find a similarly large Hn region, merging (percolating) with the first Hn region. 
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Fig. 10. — Evolution of the radiation intensity in the Lyman- Werner band, Jlw> hi units of 10 — 21 erg s _1 cm _2 str _1 Hz~ 1 . The corresponding 
photo-dissociation time is approximately ~ 100 Myrs. We consider recombination photons from Hel two-photon and Hell 

two-photon decay. 

It would be very interesting to follow the large-scale gas dynamics in this Hn complex, but the limited simulation volume 
does not allow us to study the effect of other minihalos that possibly reside further out. 

We note that the timing of the star-formation in nearby halos, and whether or not it can actually happen, is very 
sensitive to details of the configuration. In the particular simulated region, there happened to be another halo (Halo B) 
just collapsing when the first star is being turned on, at a 'favorable' location which is just outside the first Hn region. 
The survival of H2 molecules in such halos and the halo gas itself under the influence of radiation from the first star 
depends on the luminosity of the first star, distance from it, and the gas density and molecular fraction profile in the halo 
at an exactly corresponding time. Therefore, it is not straightforward to determine the fate of nearby halos in general. 
Recently, Susa & Umemura (2006) and Ahn & Shapiro (2007) addressed this issue using a number of simulations that 
explore a large parameter space. It appears that the process is indeed very complicated and is simply case-dependent. 
A 'typical' case may be found by running a large set of cosmological simulations such as those presented in the present 
paper for many halos in various environments. 



In this section, we study the thermal and chemical evolution of a gas cloud which is re-collapsing gravitationally from a 
once-ionized gas in the relic Hn region. To do so, we deliberately restrict ourselves to the vicinity of the first star; i.e., we 
ignore (possible) gas cooling and star-formation occurring in places remote from the Hn region. We also run a simulation 
without Hn region calculation, in order to make a direct comparison with runs with various radiative feedback effects 
described in the following section. 



We first consider photo-dissociation of hydrogen molecules by recombination photons from within the Hn region itself. 
The gas with ionized hydrogen and helium emits photons in a broad energy range by recombination processes. The 
dominant contribution to the Lyman- Werner band (from 11.2-13.6 eV) involves two-photon processes from helium atoms 
and ions. Since Hen is a hydrogenic ion, we can use parametric fits for the two photon profile (Nussbaumer & Schmutz 
1984); we find that each Hen 2s — > Is decay produces 0.07 LW photons. We now need to find the fraction of recombinations 
that result in two-photon decay from the 2s level, / = a e ff(2s)/ Q! totai- We can again use the fact that Hen is a hydrogen- like 
ion, so that the effective recombination coefficients obey a e g(Z,T) = Za e s(l,T/Z 2 ). Using the coefficients a c tf for HI 
tabulated in Osterbrock & Ferland (2006), we find that at T = 20, 000 K, / = 0.2, and thus that ~ 0.014 LW photons are 
produced per Hen recombination. The corresponding contribution from the 2s — > Is transition in Hel is more difficult to 
calculate since Hel is not hydrogen-like, but it is roughly comparable. 



5.2. Formation of second-generation objects in relic Hll regions 



5.2.1. Photo- dissociation by recombination photons 
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Fig. 11. — Evolution of the baryon fraction of the main halo (Halo A) and a nearby 'second-generation' halo (Halo C) from t = to t = 100 
Myrs. We define the baryon fraction as /b = Af ga s/Mtot where Mtot increases as the halo grows. 

We use directly the simulation outputs to calculate the production rate of LW photons in the relic Hn region as a 
function of time. Ideally, one would like to obtain the intensity and the spectral shape of the diffuse recombination 
radiation in the LW band at every point within the relic Hn region. This would involve a very costly radiative transfer 
calculation, whereas it is sufficient for our purposes to make a rough estimate of the importance of the recombination 
radiation in terms of molecule destruction, as is done in, e.g., Johnson & Bromm (2007). To this end, we make the 
following assumptions. We assume that the emitted LW photons are instantaneously re-distributed within a spherical 
region of radius of 1 kpc (approximately the size of the Hn region at its maximum extent). From the production rate of 
LW photons, A^lwj we approximate the intensity as 

in units of 10 _21 erg s _1 crn _2 str _1 Hz~ 1 . Note this estimate scales with the size of the assumed spherical region as 
J LW oc l/(4irR 2 ). 

Figure [TO] shows the evolution of Jlw as a function of time elapsed since the central star died. The maximum (at the 
beginning) is just about Jlw ~ 10~ 3 , and decreases slowly over several tens of million years. The large relic Hn region 
continues to expand, both by thermal pressure and by cosmological expansion, so the density in the outer part decreases, 
and helium recombination occurs progressively more slowly. 

It is clear that photodissociation by LW photons with this low flux is negligible (see equation [23] and also Machacek et 
al. 2001; Yoshida et al. 2003). As is shown in the next section, a part of the ionized gas is re-incorporated by a large mass 
halo in about 50 million years, when Jlw < 10 -5 , and thus we can safely ignore the effect of recombination radiation on 
the formation of second-generation objects in the present context. 

5.2.2. Gas evacuation and re-incorporation 

The fate of the out-flowing gas is of considerable interest, because it determines whether or not star-formation can be 
triggered promptly after the epoch of first star-formation. In other words, the gas cooling efficiency (and hence star- 
formation rate) in higher mass systems that are formed later on is critically dependent on how efficiently the halo gas is 
re-incorporated. As we show below, recollapse of the ionized gas occurs in a very complex manner. First, hierarchical 
assembly of dark matter halos is itself a complex process, and it is often found that the descendant of an early massive 
halo at one epoch is not the most massive one at a later epoch (e.g. Yoshida et al. 2003; Gao et al. 2005). We find 
exactly such a case in our simulation, where the main halo at z = 26, which is the host of the first star-forming cloud, 
grows rather slowly, and a few other smaller halos (at z = 26) grow more quickly. After the halo gas is evacuated by the 
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first star, cooling and condensation of the gas takes place at the center of another halo that is located ~ 0.5 kpc away 
from the main halo at z — 16 (see the right panel of Fig. [5]). 

In order to see clearly when this re-incorporation happens, we measure the evolution of the baryon fraction of the 
halo. The result is shown in Figure [TTJ where we normalize the baryon fraction by the global mean We show 

the evolution of baryon fraction for two halos: one is the main halo and its descendant, and the other is a different halo 
(Halo C in Fig. [9]) in which a second-generation object is formed later. Initially, the gas in the main halo is evacuated by 
radiation, with the baryon fraction rapidly decreasing from t ~ 2 Myrs when the hydrodynamic shock reaches the halo 
virial radius. The baryon fraction decreases down to a few percent(!) at ~ 40 Myrs, and starts increasing rather slowly. 
By then the halo has already grown more than twice in mass. As the gas cools, thermal pressure is lost and the gas 
eventually falls back, but full reincorporation is not achieved even after 100 Myrs. We thus conclude that the first stars, 
if they are massive, impact their surroundings by decoupling the gas distribution from that of dark matter for as long as 
a Hubble time at z ~ 20. The descendant of the main halo (Halo A) still has a very low baryon fraction of ~ 10% at that 
time. 

These results have interesting implications for the growth of putative Population III remnant black holes. Because of 
the rapid evacuation of the gas around a massive star, it is likely that the remnant black hole, even if it is formed by the 
death of the star, will not be able to accrete the surrounding gas efficiently. We defer a detailed study of this process to 
future work (Li et al., in preparation), but note briefly the physical conditions around the first star remnant. We mark a 
dark matter particle (a "black hole" ) that is located closest to the position of the central star. We then keep track of the 
position of the "black hole" and calculate the mean density and temperature with a sphere of 10 physical parsec radius 
around it. We find that the local density is ~ 0.2 cm' 3 in the beginning (see Fig. [7]), and remains at ~ 0.02 — 0.2 cm~ 3 
for over 50 million years. With this low density, the rate of gas accretion is expected to be very small. 

The strength of radiative feedback is further appreciated by the following experiment. We ran an additional simulation 
in which we include no radiative feedback from the first star. Namely, we continued the original simulation with H2 
cooling to see if star-formation can proceed if there is no feedback. The normalized baryon fraction for Halo A in this 
run remains roughly constant, at /b,norm = 0.9, and the amout of cold (T < 500K), dense (rin > 10 3 cm~ 3 ) gas increases 
from 300 M Q at z = 26 to more than 10, = OOOM by z — 18. Clearly, a large amount of cold gas should be available for 
further star- formation if there is no radiative feedback from the first star. We thus conclude that the net feedback effect 
from the first star is negative in terms of the efficiency of star-formation in the same location. 

We now return to the issue of the formation of the so-called second-generation stars. The host of the second-generation 
object, Halo C, grows quickly and so it re-captures the evaporated gas faster than Halo A does. Halo C has many 
progenitors at z = 26, which are seen as small dark matter clumps along a few filaments in the lower-right region of Halo 
A in Figure [H Merging of these numerous clumps builds up Halo C at z = 16 as is seen in the right panel of [HI The initial 
(£ < 10 Myrs) behavior of the baryon fraction for Halo C shown in Figure QT] may be misleading because there are many 
progenitors of Halo C at this epoch. Simply, the main progenitor is not well-defined then, and thus the baryon fraction 
for Halo C at t < 10 Myrs may be regarded for a halo that is one of the main progenitors of Halo C. It is also worth 
mentioning that many of the progenitors of Halo C are not within the H11 region; i.e., the gas within them is not ionized. 
Neutral and ionized gases are mixed when Halo C is assembled, although the central part of Halo C consists mostly of 
ionized gas. 

The re-captured gas starts cooling and condensing towards the center of Halo C, forming the so-called second-generation 
object. We followed the evolution of the second-generation object using the technique of Yoshida et al. (2006) until the 
central density reaches ~ 10 15 cm -3 . We emphasize that the maximum density achieved here is many orders of magnitude 
greater than previous works which study the evolution of cooling gas in Hn regions. We are, for the first time, able to 
study details of the structure and evolution of the 'cosmological' second-generation star-forming gas cloud. Figure [T2] 
shows the radial profiles around this object for density, temperature, molecular fraction (H2 and HD), and the ratio of 
enclosed gas mass to the locally estimated Bonnor-Ebert mass. 

The overall features appear similar to primordial neutral gas clouds (see, e.g., Fig. 3 of Yoshida et al. 2006). An 
important difference is the minimum temperature of ~ 50 K which is limited by Tomb = 46 K at z = 16 in the present 
calculation. The thermal evolution of the recollapsing gas is notably different from the primordial case in that H2 cooling 
is now more effective, owing to the enhanced molecular fraction. This brings the gas temperature below ~ 150 K, where 
HD cooling becomes important (see section |2~4"|) . We compare the temperature profile with that of the primordial case 
of Yoshida et al. (2006). The difference in the minimum temperature and the corresponding mass scale is clearly seen. 
Efficient HD cooling causes the temperature to bottom out at ~ 50 K (~ Tqmb) at a mass scale of 100 M®. There, 
the HD fraction increases to ~ 2 x 10 -5 , as is seen in the left-lower panel in Figure [T2l In the low temperature regime, 
the relative H2 and HD abundances approach their equilibrium values (see eauation|14|). The temperature increases 
toward the center within a mass scale of ~ IOOA/q, and the HD fraction temporarily decreases following approximately 
the equilibrium abundance (note the exponential factor exp(465K/T) in equation [2].) The density profile is close to a 
power law n oc r~ a with a ~ 2.4 depending on the distance range. The profile is slightly steeper than in the primordial 
gas clouds simulated by Abel et al. (2002), Yoshida et al. (2006) and Gao et al. (2006). Because of the lower minimum 
temperature owing to HD cooling, the temperature increase toward the center is large, by a factor of 50, and the effective 
equation of state in the collapsing gas M en ciosed < 100Mq,tiyi > 10 4 cm~ 3 , is harder than the neutral primordial case, 

2 We terminated the run at 2 = 18 because extremely small time steps are required for gas particles in cold dense gas blobs, which effectively 
stops the run. However, outputs by z = 18 are enough to make necessary comparisons. 
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P = Kp 1 with 7 ~ 1.2 — 1.3 if expressed polytropically. From the Larson-Penston similarity solution for a polytropic gas, 
the density gradient may be estimated as (Larson 1969) 

<9mp -2 

<91n7 ~ ~ 7' (29) 

which is ~ 2.5 for 7 ~ 1.2, in good agreement with the actual density gradient. The temperature increases toward the 
center within a mass scale of ~ 10M©, and the HD fraction temporarily decreases following approximately the equilibrium 
abundance. However, &o#i H 2 and HD fractions increase again when three-body reactions convert almost all the hydrogen 
atoms into molecules, making the cloud core of ~ 1M® fully molecular. The behavior of the abundances of H2 and HD 
almost parallel each other in this regime. 

To see the onset of run-away collapse, we compare the enclosed gas mass with the locally estimated Bonnor-Ebert mass 
(Bonnor 1956; Ebert 1955): 

Mbe — 



G 3/2p o l/2 ' 

« 2OM T 3 / 2 n" 1 / V V , (30) 

where mi is the first maximum mass of the solution for the isothermal Lane-Emden equation (see, e.g. Stahler & Palla 
2004), and \i and 7 denote the mean molecular weight and adiabatic index, respectively. We approximate the external 
pressure by its local value taken from the radial density and temperature profiles. (Note that Mbe evaluated in this 
manner is essentially the same as the local Jeans mass.) 

We find that the enclosed gas mass exceeds the Bonnor-Ebert mass at a mass scale of M ~ 40M Q . This is the 
characteristic mass of the collapsing gas cloud. The so-called second-generation object formed in this process may be a 
formation site of low-mass primordial stars (Mackey, Bromm & Hernquist 2003; Johnson & Bromm 2006). The mass of 
such stars and the possible difference from the first stars are of particular importance. Until the time of the last output, 
when the central density is ~ 10 15 cm~ 3 , we observe no sign of fragmentation of the pre-stellar cloud, and thus it is likely 
that a single proto-stellar seed is formed in this object. Detailed calculations of proto-stellar evolution for this object will 
be presented elsewhere (Yoshida, Omukai & Hernquist, in prepartion). 
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Fig. 12. — Radial profiles around the 'second-generation star' for density, temperature, molecular fraction, and the ratio of enclosed gas 
mass to locally estimated Bonnor-Ebert mass. The density is plotted as a function of distance from the center, whereas the other three 
quantities are plotted as a function of enclosed gas mass. The dashed line in the density profile plot indicates a power law of oc r~ 2,4 . The 
horizontal dashed line in the temperature profile plot indicates the CMB temperature at z = 16. We also show the temperature profile of the 
primordial proto-star in Yoshida et al. (2006) for comparison (dot-dashed line). In the lower-left panel, we shift the HD fraction two decades 
upward to simplify comparison with H2 fraction. 
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Fig. 13. — Mass-weighted mean H2 fraction (left) and the normalized baryon fraction (right) against virial temperature for halos in the 
relic Hn region at z=16. Filled circles are halos in the relic Hn region in the run without background radiation, whereas open circles are those 
in the run with Jlw=0-1. We also show the molecular fractions and baryon fractions in the run without Hn region calculation (triangles). 
The dashed line in the top panel is the asymptotic H2 fraction in a collapsed neutral gas oc T 1,52 derived by Tegmark et al. (1997). The big 
square in the bottom panel is Halo C, which has a mixture of (once-)ionized gas and neutral gas. 



5.3. The effect of ultra-violet radiation background 

So far, we have not considered external radiation from other stars, assuming that the first star in the simulation is the 
very first radiation source. We have also shown that recombination radiation from the relic Hn region itself is unimportant. 
In a more general cosmological context, there may be an early UV background radiation built up by stars/galaxies that are 
formed elsewhere. In particular, the intergalactic medium is nearly optically thin to photo-dissociating far-UV radiatiorQ, 
and thus it is conceivable that an early far-UV background radiation is quickly built up when stars/galaxies are formed. 
It is still unclear whether the net sign of feedback is positive (Ricotti, Gnedin & Shull 2002) or negative (Oh & Haiman 
2003) in relic Hn regions relative to primordial neutral regions. The strength of the background LW radiation appears to 
play a key role in modulating between these regimes. Mesinger, Bryan & Haiman (2006) recently examined these issues, 
and found that there is a critical value of Jlw ~ 0.01 x 10 _21 ergs _1 cm _2 str _1 Hz _ for which H2 cooling in relic Hn 
regions is strongly suppressed. However, their simulations do not take gas self-shielding into account, and thus the critical 
value they obtained is likely to be an under-estimate. 

In principle, the background radiation itself evolves and thus the intensity, as well as the spectrum, changes as a function 
of time. Since the relevant energy range for photo-dissociation is narrow, between 11.18 eV and 13.6 eV0, we assume 
that the spectrum is flat in this range, and we set the intensity to be J(v) = 0.1 x 10 _21 ergs~ 1 cm _2 str _1 Hz _1 . We also 
implement the effects of H2 self-shielding, as described in section 13.51 

We evolve the system down to z — 16 with the radiation on and compare the outputs from those without background 
radiation. Figure [TBI shows the mean molecular fraction in halos that are in the relic Hn region. In the figure, open circles 
are for halos in the run with a LW radiation field, and solid circles are for halos in the other run. For reference, we also 
show the molecular fractions in the run without the Hn region calculation described in the previous section (triangles). 
The dashed line in the top panel in Figure [T3] is the asymptotic molecular fraction in a neutral primordial gas as calculated 
in Tegmark et al. (1997), which scales as oc T . The molecular fraction in the run without radiative transfer calculation 
is well described by the model. The so-called positive feedback in relic Hn regions is clearly seen as enhanced large 
molecular fractions in the case without LW radiation. However, the positive feedback is effective only if the background 
far-UV radiation is weak. Even though molecule formation is promoted in relic Hn regions, UV background radiation 
with a modest intensity is enough to prevent the gas from cooling and collapsing by molecular hydrogen cooling. This 
is explained by the low-density, and equivalently high entropy, of the gas in the relic Hil region. Oh & Haiman (2003) 
argue that photo- ionization heating increases the gas entropy, which prevents gas condensation in low-mass halos, possibly 
acting as negative feedback. 

3 The primordial gas at very high redshift is not optically thin, because there is a trace amount of H2 molecules left over in the post- 
recombination epoch (e.g., Ricotti et al. 2001). However, the residual H2 molecules in the diffuse intergalactic medium are quickly destroyed 
by the first luminous sources. 

4 Photo-detachment of H — by lower energy photons than LW photons could also affect the formation of H2 molecules, and hence it should, 
in principle, be taken into account. However, the relative intensities between the relevant energy ranges are not well-constrained because of 
the uncertainties in the sources of background radiation. We thus ignore radiation below the LW bands as well, in order to isolate the effect of 
photo-dissociation in the present paper. 
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Fig. 14. — Radial profiles around the 'second-generation star' in the buried Hn region for density, temperature and molecular fraction. The 
lower-right panel shows the gas distribution within the halo (the side length is J? v ir = 80pc) and a close-up on the central 0.4 pc region. The 
dashed line in the density profile plot indicates a power law of oc r~ 2,3 . 

The right panel of Figure [T3l shows the normalized baryon fraction for the 'ionized' halos in the relic Hn region. Again, 
we show the baryon fractions for halos in the run without the Hn region calculation (triangles). Small-mass halos in the 
relic Hn region have a significantly lower baryon fraction, and hence the mean gas density within them is very small. 
Consequently, molecule formation does not proceed efficiently, and self-shielding does not counter-act photo-dissociation 
for the adopted intensity of the background radiation. Halos in the relic Hn region have a significantly lower baryon 
fraction, and hence the mean gas density within them is very small. Consequently, molecule formation does not proceed 
efficiently, and self-shielding does not counter-act photo-dissociation. Oh & Haiman (2003) estimate that the gas fraction 
in mini-halos should be as small as ~ O.f , which is in good agreement with the results of our three-dimensional simulations 
with hydrodynamics, as is seen in Figure fTSI The mean molecular fractions are below 10~ 5 in the case with radiation, much 
less than the critical fraction necessary for gas cloud collapse (Yoshida et al. 2003; Oh & Haiman 2003). It is interesting 
that the fractions are also smaller than the asymptotic fraction (dashed line) for a primordial gas, which indicates a 
net negative feedback effect in comparison with early primordial gas clouds. We have set the radiation intensity to be a 
constant. For a weaker UV radiation field, the result will be close to the no radiation case, whereas for a stronger radiation 
field, photo-dissociation is more effective and thus the molecular fraction will likely be smaller than shown in Fig. [T31 
Note, however, that the formation and destruction processes are highly non-linear, and hence the resulting molecular 
fraction will not simply scale with the radiation intensity. 

We close this section by mentioning that the effect of meta-galactic ionizing radiation will be more profound in regulating 
star- formation in galaxies at high redshift. A number of authors addressed this issue using one-dimensional (Thoul & 
Weinberg 1996; Kitayama & Ikeuchi 2000) and three-dimensional (Susa & Umemura 2004a, b; 2006) calculations. These 
authors show that dwarf galaxy formation is regulated by both global UV background radiation and local UV sources. 
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Once global reionization is completed, formation of H2 molecules is substantially suppressed and thus it is likely that 
star-formation takes place only in high-mass halos in which the gas cools by atomic processes. As we have shown in 
this section, photo-dissociating radiation with a modest intensity is enough to delay star-formation in relic Hn regions 
at least for over tens of million years. If ionizing photons by global or local luminous UV sources arrive soon after the 
photo-dissociating photons, star-formation in the region will be further delayed. 

5.4. Buried HII regions 

We have studied the formation and evolution of an Hn region that extends much further than the host halo's virial 
radius. There may be a distinct case that the I-front is trapped well inside the host halo, making a compact Hn region. 
Namely, if the I-front remains D-type during the stellar lifetime, the size of the Hn region would be much smaller than 
the virial radius. This corresponds to a case where the central star is less luminous, or the host halo mass is very large. 
KYSU simulated a case with a high mass halo, Mh a i = 1O 7 M , and found that a D-type front travels only up to ~ 2pc, 
and then falls back to the center. It is interesting to study the evolution of such 'buried' Hn regions and the formation of 
second-generation stars in them. 

To model this case, we use the same simulation as in the previous sections but switch off the star when the I-front 
reaches ~ 5 pc from the center. By this time, the central ~ 4000M© part has been ionized. The ionized gas still has a 
high density and thus recombines quickly, while expanding because of its thermal pressure, and the central part evolves 
essentially as in the isobaric case presented in section 2. The mean gas density at the center is as large as ~ 3 cm" , 
and the thermal energy is lost in less than a few million years mainly by hydrogen line cooling. The gas re-collapses and 
condenses in about ten million years, and then a dense cold gas blob is formed at the center. We again use the technique 
of Yoshida et al. (2006) to follow the gas evolution to densities of ~ 10 15 cm -3 . Figure [T4l shows the density, temperature 
and molecular fraction profiles around the gas cloud. The lower-right panel shows the gas distribution within a volume of 
100 pc on a side, and within the central 0.4 pc region. The extent of the compact Hn region can be seen as the central, 
roughly spherical part colored in red. The radial profiles look similar to the second-generation object discussed in the 
previous section, except that the minimum temperature is slightly higher, around 100K. This is because the collapsing 
cloud has a mixture of an ionized gas and a neutral gas. As soon as the central radiation source is switched off, the sharp 
boundary (a pressure jump) of the buried Hn region drives significant mixing of the gas in both the ionized and neutral 
sides. The two gases collapse toward the potential well of the host dark halo together, and thus the final collapsing gas 
cloud contains a substantial amount of neutral gas in and around the envelope. HD cooling does not operate in such a 
neutral gas packet, maintaining its temperature just above 200K. When we calculate average temperatures over spherical 
shells, we obtain a slightly higher minimum temperature because of the 'warmer' neutral gas packets. By calculating 
the ratio of the enclosed gas mass to the locally estimated Bonnor-Ebert mass (equation [3D])j we find that the mass of 
the cloud that is run-away collapsing is ~ 80M Q , roughly consistent with the cloud mass in the re-collapse calculation 
in section 15.2.21 Throughout the evolution, we observe no sign of fragmentation of the cloud up to the final output time 
when the cloud core density reaches ~ 10 15 cm -3 . Following the argument of Omukai & Palla (2003) and Yoshida et al. 
(2006), we expect the central proto-star in this buried Hn region will grow, to be a massive metal- free star. 

6. SUMMARY AND DISCUSSION 

We have studied the formation and structure of early cosmological Hn regions using three-dimensional radiation hy- 
drodynamics simulations. We have also examined in detail the formation of second-generation stars, by following cooling 
and condensation of primordial gas in relic Hn regions. With the three-dimensional simulations, we can eliminate a 
number of uncertain assumptions regarding the formation of early Hn regions. First, the simulations allow us to directly 
calculate the ionizing photon escape fraction from cosmological minihalos. With enough mass resolution, our calculation 
self-consistently includes the effect of gas clumping. We obtain photon escape fractions of order unity, for the cases with 
rather massive (M* > 60M Q ) central stars. While our radiation tranfer scheme is based on a few approximations, the 
accuracy of the method is tested against the results of previous one-dimensional calculations (KYSU; Whalen et al. 2004). 
We conclude that the technique is well-suited for the particular problem of early Hn regions. 

We find that a large volume of a few kiloparsec diameter is ionized by a single massive star, within which a smaller 
Hem region is embedded. The ionized hot gas is evacuated at a velocity of ~ 30 km/sec. After the central star dies off, 
the gas recombines, and cools first by atomic line cooling, then by H2 cooling, and finally by HD line cooling down to a 
few tens Kelvin. The fractionation of HD/H2 occurs when the gas temperature becomes as low as 100 K, and the HD 
fraction reaches 10~ 5 . Run-away collapse is triggered when the cloud mass exceeds ~ 40M Q at the temperature minimum 
of ~ 40 — 50K. It is substantially smaller than the corresponding mass scale for the first generation stars, ~ 300M Q 
(Yoshida et al. 2006), indicating that second-generation stars formed under these conditions likely have smaller masses 
than the first stars. At least the maximum mass will be limited by the gas clump mass. It is worth noting that the 
elemental abundance patterns of hyper metal-poor stars recently discovered by Christlieb et al. (2002) and Frebel et al. 
(2005) indicate that the early metal-enrichment was caused by supernova explosions with a progenitor mass of ~ 25Mq 
(Umeda & Nomoto 2003; Iwamoto et al. 2005). We speculate that the second-generation stars formed in relic Hn regions 
likely have a similar mass. 

It is important to note that the re- incorporation and recollapse of the ionized gas happens very late. Full re- incorporation 
of the gas is not achieved even after a hundred million years, but a small amount of the gas at the center of a growing 
halo recollapses to form a star-forming gas cloud in a hundred million years. O'Shea et al. (2005) report a shorter time 
scale for recollapse. The discrepancy can be understood from the fact that, in their calculation, the density structure is 
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unchanged during the evolution of the Hn region. The highest density regions have the shortest recombination times and 
radiative cooling times, and thus the gas evolves faster than it should if hydrodynamics and radiation are coupled. By 
employing a more accurate treatment of radiation transport, we are able to obtain a realistic estimate for the timescale 
of suppression of star-formation in early halos. We also examined the feedback effect of external UV radiation, and found 
that it strongly reduces the H2 fraction and baryon fraction in low-mass halos. Our results suggest that the long-term 
evolution and the fate of the ionized gas in and outside relic Hn regions are determined by a rather complex combination 
of many factors such as the source luminosity, relative positions, exact timing of collapse, strength of background radiation 
etc., and hence it is very difficult to establish a 'rule' of the sign of the net feedback effect. This is also supported by 
recent numerical studies by Ahn & Shapiro (2007) and Susa & Umemura (2006) who explored a large parameter space 
in order to examine the net effect of photo-dissociation. We argue that it is necessary to carry out a large cosmological 
simulation with all these relevant physics implemented, as done in the present paper, in order to to follow the evolution 
in relic Hn regions; i.e., probable sites of (proto-)galaxy formation. 

In the present paper, we did not consider cooling by heavy elements, assuming that the gas in Hn regions remains 
chemically pristine although its ionization fraction changes. Metal-enrichment can be caused not only by supernova 
explosions but also by stellar winds during the evolution of the central massive star via mixing and self-enrichment (e.g., 
Vink & de Koter 2005). It is thus conceivable that at least a small amount of carbon can be dispersed in the central part 
of the Hn region. It is well known that trace amounts of heavy elements such as carbon and oxygen can significantly 
change the gas cooling efficiency (e.g., Bromm & Loeb 2004; Omukai ct al. 2005; Santoro & Shull 2006). Intriguingly, 
however, Omukai et al.'s calculations show that the thermal evolution of a collapsing gas with a modest metal content, 
but without dust, is similar to that found in our simulation of an ionized primordial gas with HD cooling. Although 
radiative cooling by the Cn fine structure transition can bring the gas temperature as low as 10 K, whereas HD cooling 
is efficient at T > 30 K, the minimum gas temperature is limited by Tomb > 30 K at z > 10, and thus gas-phase metal 
enrichment may not substantially affect the results of our calculations. We argue that the existence of dust may play a 
more important role in significantly changing the characteristic mass of collapsing gas clouds. The gas temperature can 
be kept low at high densities by dust thermal emission, which brings the collapse mass scale to the 'opacity-limited' value 
of < IM Q (Low & Lynden-Bell 1976; Schneider et al. 2006). 

We now discuss the implications of our results for early galaxy formation. While the very first generation stars are 
likely to be born in low-mass dark halos, it is often argued that efficient star formation takes place only in large halos 
in which the gas cools via atomic processes (e.g., Barkana & Loeb 2000). A crude assumption in such analyses is that 
the gas distribution traces that of dark matter; i.e., the baryon fraction is close to the cosmic mean, so that the gas 
cools efficiently at the centers of dark halos. Our simulation results invalidate these assumptions. If the first stars are 
massive (not necessarily very massive) , they decouple the gas distribution from that of dark matter for as long as a Hubble 
time at z ~ 20. We argue that efficient star-formation does not take place immediately in dwarf (proto-)galaxies whose 
progenitors have hosted massive PopIII stars earlier. 

The first stars may leave remnant black holes at the end of their lives, unless they trigger complete disruption by pair- 
instability supernova explosions. Efficient accretion onto the remnant black holes could lead to formation of miniquasars 
(Kuhlen & Madau 2005). Our simulation results indicate, however, that accretion is rather inefficient at least for some 
tens of millions years because the progenitors, massive stars, emit a large number of UV photons and the surrounding 
gas is evacuated out of the gravitational potential well of the host halo. It will be interesting to simulate the inefficient 
accretion process onto the remnant Population III black holes. 

The first Hil regions may imprint a distinct signature as a source of CMB secondary anisotropics via the kinematic 
Sunyaev-Zel'dovich effect. The diameter of a few kiloparsec proper region at z = 10 — 20 is about a few arcsecs in angular 
size. Hence, if a large number of early Hn regions are formed, they will produce CMB temperature fluctuations, at arcsec 
scales in the angular power spectrum. If the first Hn regions are clustered in cosmological high-density peaks such as 
those studied by Reed et al. (2005), they will make a large Hn bubble with a mega-parsec diameter at z > 10, which may 
possibly be detected by ALMA with multiple channels with a very long-time exposure. Finally, it has been long known 
that neutral hydrogen at high redshift may be directly detectable through the redshifted 21cm line (see Furlanetto, Oh 
& Briggs (2006) for a recent review). This can be seen either in absorption or emission against the CMB. We have made 
a brightness temperature map following Nusser (2005) using the simulation outputs. We find that the overall signal is a 
strong function of time (equivalently the evolutionary phase of the Hn region) . While the Hn region during the central 
stellar lifetime appears essentially as a "void" in the map, because the HI content is extremely small, the relic Hn region 
can be seen as a strong emission source with the differential temperature against the CMB up to 6Tb ~ +50mK for an 
angular resolution of ~ 0.02". While it is unlikely that currently planned radio telescopes can directly detect individual 
Hn patches, the fluctuations in the 21cm background from early Hn regions can be significantly enhanced just like those 
from cosmological minihalos (e.g., Shapiro et al. 2006). We will study the observability of early Hn regions in future 
work. 
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APPENDIX 

Chemical reactions involving Hjj~ and HcH + that may be important in a primordial gas are summarized below. 
For H^", the main formation paths are 

H++H 2 ^H++H, (1) 

and 

H+ + H^H++ 7 . (2) 
We include the reverse reaction of these and dissociative recombination 

H+ + e -> H 2 + H, (3) 
- 3H. (4) 

For HeH + molecules, the main formation paths are the radiative association 

Hc + H 2 -> HeH+ +H+, (5) 

and the inverse rotational predissociation 

He + H + -► HcH+ + 7. (6) 
We include the reverse processes and the dissociative recombination 

HeH+ + e -> He + H, (7) 

and a formation path of H^ 

HcH + + H 2 -> H3 + He (8) 

The reaction rates are taken from Stancil et al. (1998). Hirata & Padmanabhan (2006) recently revised post-recombination 
calculations of the early Universe. They argue that H^ formation via 

HcH+ + H -> H+ + He (9) 

is the dominant path at high redshifts (z > 200) . We have checked that this process does not affect the Hj abundance in all 
the thermal phases and evolution considered in this paper. Also, photo-dissociation by the cosmic microwave background 
radiation is unimportant in the redshift range we consider. 

We have run the same calculations of an isobarically cooling gas as in section 12.41 including the above species and 
reactions. We conclude that the number fractions of these ionic molecules are extremely small, and do not affect the 
thermal evolution of the gas. Cooling by H^ ions can potentially be important in collapsing gas clouds at very high 
densities, njj ~ 10 7 — 10 10 cm~ 3 when there is some (weak) ionization source (S. Glover, private communication), but it 
is outside the regime and conditions we consider in the present paper. 
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Table 1 

REACTION RATE COEFFICIENTS: Deuterium Chemistry 



Reactions Rate Coefficients (cm 3 s ) Reference 

(Dl) T>+ + e^D + hv fc D i = 3.6 x 10- 12 (T/300)-° 75 1 

(D2) D + H+^D++H fc D2 = 2.0 x 10- 10 T a402 exp(-37.1/T) -3.31 x 10- 17 T 148 2 

(D3) D++H^D + H+ fc D3 = 2.06 x 10- 10 T - 396 cxp(-33.0/T) + 2.03 x I0- 9 T-°- 332 2 

(D4) D + U^UD + hv fc D4 = 1.0 x 10~ 25 3 

(D5) D + H 2 ^H + HD fc D5 = 9.0 x 10" 11 cxp (-3876/T) 4 

(D6) HD++H^H++HD fc D6 = 6.4 x lO" 10 3 

(D7) D++H 2 ^H++HD fc D7 = 1.6 x lO" 9 4 

(D8) HD + H^H 2 +D fc D8 = 3.2 x 10" 11 cxp (-3624/T) 3 

(D9) HD + H+^H 2 + D+ fc D9 = 1.0 x 10~ 9 cxp (-464/T) 3 

(40) D + H+ -> HD+ + hv fc D io = dex[-19.38 - 1.523 logT+ 1.118(logT) 2 - 0.1269(logT) 3 ] 3 

(Dll) D+ + H -> HD+ + hv fc D ii = dcx[-19.38 - 1.523 logT+ 1.118(logT) 2 - 0.1269(logT) 3 ] 3 

(D12) HD++e^H + D k D12 = 7.2 x lO" 8 ^ 1 / 2 3 

(D13) D + e -> D" + hv k D13 = 3.0 x 10- 16 (T/300) 95 cxp(-T/9320) 1 

(D14) D+ + D" -> 2D fcou = 5.7 x 10- 8 (T/300)-°' 5 1 

(D15) H++D-^D + H fc D15 = 4.6 x 10- 8 (T/300)-°- 5 1 

(D16) H-+D^H + D- fcoie = 6.4 x 10- 9 (T/300) - 41 1 

(D17) D-+H^D + H- fc D17 = 6.4 x 10- 9 (T/300) ' 41 1 

(D18) D"+H^HD + e fc D18 = 1.5 x 10- 9 (T/300)- 1 1 



Note.— (1): Galli & Palla (1998); (2): Savin (2002); (3): Stancil, Lepp & Dalgarno (1998) (4): Wang & Standi 
(2002); 



